home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 / Ham Radio 2000.iso / ham2000 / antenna / yagiu112 / solve.c < prev    next >
C/C++ Source or Header  |  1995-08-21  |  5KB  |  130 lines
  1. #include <stdio.h>
  2. #include <memory.h>
  3. #include <errno.h>
  4. #include "yagi.h"
  5. extern int errno;
  6.  
  7. #define MAX_NON_ITERATIVE 100000
  8.  
  9. void solve_equations(double frequency, int driven, int parasitic, double **driven_data, double **parasitic_data, double *v, double **z, double *pin, struct FCOMPLEX *voltage, struct FCOMPLEX *current, struct FCOMPLEX *input_impedance, struct element_data *coordinates, double **A, double *b, int *indx)
  10. {
  11.     int elements, element_number, i, j;
  12.     double d, **A_pre_LUdcmp, *b_copy, t1, t2;
  13.     struct performance_data rubbish;
  14.     struct flags flag;
  15.     elements=driven+parasitic;
  16.     fill_z_matrix(frequency,driven,parasitic,driven_data,parasitic_data, z);
  17.     /* Now we must fill the V vector, which gives the voltage at the centre of
  18.     each driven element. Since we know the magnitude of the voltage and phase,
  19.     we can calculate the real and imaginary components of voltage. NB */
  20.     fill_v_vector(driven, parasitic, driven_data, v);
  21.  
  22.     /* We now have the voltage vector V, and the impedance matrix Z. All
  23.     we need to do now is solve a set of NxN equations, where N is the
  24.     number of elements, to get N values of current, which we can put in
  25.     the I matrix.
  26.     Unfortunately, this is not trivual, by any standards!!!  It is
  27.     complicated by the fact the the N equations are all complex. 
  28.     The method used us a brute-force approach mentioned by Press
  29.     in their 2nd Edition of the Numerical Recipes in C book. Here we
  30.     put the data in the form:
  31.     
  32.     |Zr -Zi|   |Ir|   = |Vr|   which we will call A.x=b
  33.     |Zi  Zr|   |Ii|     |Vi| 
  34.         
  35.     so the Z data goes now in a 2Nx2N matrix. This is a bit wasteful of space
  36.     and time, but it will do here. */
  37.     /* Copy impedance data from 'z' to matrix A */
  38.     copy_complex_data_to_real_matrix(elements,z,A);
  39.     /* The following function prints to stantard output the z matrix of the
  40.     antenna. It is really only used during debugging, so it can be commented
  41.     out normally. */
  42.     /* print_z_matrix(frequency,elements,z);  */
  43.     /* read voltage data from v into b */
  44.     for(i=1;i<=elements;++i)
  45.     {
  46.         b[i]= v[2*i-1];           /* real data */
  47.         b[i+elements]=v[2*i];    /* imaginary data */
  48.         voltage[i].r=v[2*i-1];
  49.         voltage[i].i=v[2*i];
  50.     }
  51.     /* If there are a lot of elements, its possible that rounding errors
  52.     will destroy the accuracy of the results. If there are more than 
  53.     MAX_NON_ITERATIVE elements, we do the usual LU decompositioon and
  54.     back-substitution, then do an iterative improvement of the 
  55.     solution */
  56.     if(elements > MAX_NON_ITERATIVE)
  57.     {
  58.         A_pre_LUdcmp=dmatrix(1L, 2L*(long)elements, 1L, 2L*(long) elements);
  59.         b_copy=dvector(1L, 2L*(long)elements);
  60.         /* I had troubled using memcpy to copy matrices - probably becuase
  61.         they are offset by  1 */
  62.         copy_matrix(2*elements, 2*elements, A_pre_LUdcmp, A);
  63.         /* copy vector b */
  64.         for(j=1;j<=2*elements;++j)
  65.             b_copy[j]=b[j];
  66.     }
  67.     /* Perform a LU decompositon of A */
  68.     ludcmp(A, elements*2, indx, &d);
  69.     /* We now have the voltages in b. After lubksb is run, we get the 
  70.     currents in b . A contains an LU decomposed version of the original A*/
  71.     lubksb(A, 2*elements, indx, b); /* current's in b after lubksb has run*/
  72.     if(elements>MAX_NON_ITERATIVE)
  73.     {
  74.         t1=b[1];
  75.         mprove(A_pre_LUdcmp, A, 2*elements, indx, b_copy, v);
  76.         free_dmatrix(A_pre_LUdcmp,1L, 2L*(long)elements, 1L, 2L*(long) elements);
  77.         /* free_dvector(bcopy,1L,2L*elements); */
  78.        t2=b[1];
  79.         if(t1!=t2)
  80.         {
  81.             printf("b[1]'s differ before and after mprove: before = %.16lf after =%.16lf\n", t1, t2);
  82.             exit(1);
  83.         }
  84.     }
  85.     /* Put currents an FCOMPLEX matrix current[element] 
  86.     currents are stored in b as r,r,r,r ..... i,i,i,i */
  87.     for(element_number=1;element_number<=elements;element_number++) 
  88.     {
  89.         current[element_number].r=b[element_number];
  90.         current[element_number].i=b[element_number+elements];
  91.     } 
  92.     z_input(voltage[1], current[1], input_impedance);
  93.     *pin=calculate_power_input(input_impedance->r,current[1]); /* in Watts */
  94.     /* Sometimes the input power goes < 0. I'm not sure why this occurs, but
  95.     the implication is that the antenna is a source of energy, which is silly.
  96.     The gain is negative then (not in dB), so taking the log to put in dB
  97.     would cause a numerical error. When this occurs, we print a message to
  98.     stderr, then write the data in a file "problem". The user may then 
  99.     investigate what went wrong. */
  100.     if(*pin <= 0.0)
  101.     {
  102.         /* fprintf(stderr,"Input power less than 0! Undesirable, but non-fatal error.\n"); */
  103.         *pin=1e20; /* Force it to a large value */
  104.    /* do_since_better(0,"problem","problem",*input_impedance, \
  105.     rubbish,flag, "This causes Pin < 0.0 - findout why!",frequency, frequency, \
  106.     frequency,frequency/10.0, elements, driven, parasitic, \
  107.     180.0, driven_data, parasitic_data, 1.0, 0.0); */
  108.     } 
  109.     for(i=1;i<=driven;++i)
  110.     {
  111.         coordinates[i].x=driven_data[i][X];
  112.         coordinates[i].y=driven_data[i][Y];
  113.         coordinates[i].length=driven_data[i][LENGTH];
  114.     }
  115.     for(i=driven+1;i<=elements;++i)
  116.     {
  117.         coordinates[i].x=parasitic_data[i-driven][X];
  118.         coordinates[i].y=parasitic_data[i-driven][Y];
  119.         coordinates[i].length=parasitic_data[i-driven][LENGTH];
  120.     }
  121.  
  122. #ifdef DEBUG
  123.     if(errno)
  124.     {
  125.     fprintf(stderr,"Errno =%d in solve.c\n", errno);
  126.     exit(1);
  127.     }
  128. #endif
  129. }
  130.